Escaping free-energy minima 
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We introduce a novel and powerful method for exploring the properties of the multidimensional 
free energy surfaces of complex many-body systems by means of a coarse-grained non-Markovian 
dynamics in the space defined by a few collective coordinates. A characteristic feature of this 
dynamics is the presence of a history-dependent potential term that, in time, fills the minima in the 
free energy surface, allowing the efficient exploration and accurate determination of the free energy 
surface as a function of the collective coordinates. We demonstrate the usefulness of this approach 
in the case of the dissociation of a NaCl molecule in water and in the study of the conformational 
changes of a dialanine in solution. 
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Molecular dynamics and the Monte Carlo simulation 
method have had a very deep influence on the most di- 
verse fields, from materials science to biology, and from 
astrophysics to pharmacology. Yet in spite of their suc- 
cess, these simulation methods suffer from limitations 
that reduce the scope of their applications. A severe 
constraint is the limited time scale that present-day com- 
puter technology and sampling algorithms explore. In 
particular, there are many circumstances where the free 
energy surface (FES) has several local minima separated 
by large barriers. Examples of these situations include 
conformational changes in solution, protein folding, first 
order phase transitions, and chemical reactions. In such 
circumstances a simulation started in one minimum will 
be able to move spontaneously to the next minimum only 
under very favourable circumstances. A host of methods 



have been suggcs 
explore the FESi 
ize the transition stat 




restriction and to 
or to character- 



Here we find a new solution 
to this pr pbJep i by combi ping , the ideas of coarse-grained 
dynamicsEJIiZi on thp-FESliErES with those of adaptive bias 
potential methodsHU obtaining a procedure that allows 
the system to escape from local minima in the FES and, 
at the same time, permits a quantitative determination 
of the FES as a by-product of the integrated process. 



I. METHODOLOGY 

We shall assume here that there exist a finite num- 
ber of relevant collective coordinates s^, i — l,n where 
n is a small number, and we consider the dependence of 
the free energy T{s) on these parameters. Practical ex- 
amples of appropriate choices of these variables will be 
given below. The exploration of the FES is guided by 
the forces Ff = —d!F/ds\. In order to estimate these 
forces efficiently, we introduce an ensemble of P repli- 
cas of the system, each obeying the constraint that the 
collective coordinates have a pre-assigned value Si — s\, 
and each evolved independently at the same temperature 
T . Since the P replicas are statistically independent, the 
estimate of thermodynamic observables (e.g. the forces 
on the constraints) is improved with respect to an eval- 



uation on a single replica, and it can be parallelized in 
a straightforward manner. The constraints are imposed 
on each replica via Ijlie standard methods of constrained 
molecular dynamicaia by adding to the Lagrangean a 
term X]i=i n '^'(■^^ ~ ^\)^ where \i are Lagrange multi- 
pliers. Averaging over the time and over the replicas we 
can evaluate the deriyative of the free energy relative to 
the s*'s as — {Xi)p'a. For the sake of simplicity we ne- 
glect here the small kinematic correction term discussed 
in Ref.Q, although this can easily be added. Taking a cue 
from the work of Gear and Kevrikidialj, we use these 
forces, determined by a microscopic dynamics, to define 
a coarse-grained dynamics in the space of the s^'s. In 
our case the definition of this dynamics is rather arbi- 
trary and designed only to explore the FES efficiently. 
The dynamics is defined from the discretized evolution 
equation: 
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In eqn. (^ we have introduced the scaled variables 
(7* = sj/Asi and the scaled forces (pj = FfAsi, where 
Asi is the estimated size of the FES well in the direction 
Si, \(f>^\ is the modulus of the n-th dimensional vector (/)* 
and Sa is a dimensionless stepping parameter. It should 
be stressed that eqn. (^ has the form of a steepest de- 
scent equation in the direction given by the forces (pi, and 
does not imply any real dynamical evolution. The index 
t is only used to label the states. After the collective co- 
ordinates are updated using eqn. (|l|), a new ensemble of 
replicas of the system with values cr,*"''^ is prepared, and 
new forces F^'^^ are calculated for the next iteration. At 
the same time the driving forces are evaluated from the 
microscopic Hamiltonian in short standard microscopic 
molecular dynamics runs. Since there is no dynamical 
continuity between the P replicas at different iterations, 
one can use large values of Sa and move very efficiently 
in the space of the collective coordinates. 

Clearly eqn. ([^) alone cannot guarantee an efficient 
exploration of the FES, nor is it useful to determine the 
FES. This task can be achieved if we repJ*,G^ the forces 
in eqn. (nT) with a history-dependent terr 
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where the height and the width of the Gaussian W and 5a 
are chosen in order to provide a reasonable compromise 
between accuracy and efficiency in exploring the FES, as 
we will show below. The component of the forces coming 
from the Gaussian will discourage the system from revis- 
iting the same spot and encourage an efficient exploration 
of the FES. As the system diffuses though the FES, the 
Gaussian potentials accumulate and fill the FES well, al- 
lowing the system to migrate from well to well. After a 
while the sum of the Gaussian terms will almost exactly 
compensate the underlying FES well. 




FIG. 1: Time evolution of the sum of a one-dimensional 
modelpotential V((t) and the accumulating Gaussian terms of 
eqn. (bl) . The dynamic evolution (thin lines) is labelled by the 
number of dynamical iterations (^). The starting potential 
(thick line) has three minima and the dynamics is initiated in 
the second minimum. 

A typical example of this behaviour is shown in Figure 
1, in which a dynamics in eqn. (|l|) is used to explore 
a one-dimensional potential energy surface (PES) V{(t) 
with three wells. The dynamics starts from a local mini- 
mum that is filled by the Gaussians in 20 steps. Then 
the dynamics escapes from the well from the lowest en- 
ergy saddle point, filling the second well in ~ 80 steps. 
The second highest saddle point is reached in ~ 160 steps, 
and the full PES is filled in a total of ~ 320 steps. Hence, 
in the case of this example, since the form of the poten- 
tial is known, it can be verified that, for large t and if the 
width of the Gaussians is sufficiently small with respect 
to the length of a typical variation of V, 
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modulus an additive constant. The same property can 
be verified for any multidimensional potential if the vari- 
ables <Ti are allowed to vary in a finite region. Hence, 
we assume that eqn. (H) holds also if the method is 
used to explore the FES, and that the free energy can 
be estimated from eqn. (|^) for large t. This procedure 
resembLea the algorithm recently proposed by Wang and 
LandaulliJ, in which a time-dependent bias is introduced 
in order to modify the density of states to produce locally 
fiat hystograms. 

We observe an interesting isomorphism between our 
dynamics and the ordered deposition of the beads of a 
polymer chain on the FES. In fact we can regard cr* as the 
position of a monomer in the n-dimensional parameter 
space. Each monomer is held at a distance 5a from the 
previous one and monomers repell each other with the 
Gaussian term of eqn. (H). Neglecting terms of order 
(5ct^ each monomer is placed by eqn. in the position 
of free enegy minimum compatible with the requirement 
of minimum overlap with previous beads. Hence, the 
polymer chain generated in this manner fills the wells 
in the FES. This isomorphism helps to clarify how our 
approach works and why it can be made more precise by 
reducing 5a. 

The efficiency of the method in filling a well in the 
PES (or in the FES) can be estimated by the number of 
Gaussians that are needed to fill the well. This number is 
proportional to (l/5cr)", where n is the dimensionality of 
the problem. Hence, the efficiency of the method scales 
exponentially with the number of dimensions involved. If 
n is large, the only way to obtain reasonable efficiencies 
is to use Gaussians with a size comparable to that of the 
well. On the other hand, a sum of Gaussians can only 
reproduce features of the FES on a scale larger than ~ 5a. 
A judicious choice of Asi, W and 5a will ensure the right 
compromise between accuracy and sampling efficiency, 
and the optimal height and width of the Gaussians are 
determined by the typical variations of the FES. 

If prior information on the nature of the free-energy 
well is not available, the scaling parameters As, are cho- 
sen by performing short coarse-grained dynamic runs 
without bias potential (see eqn. (^). In such a case 
the system moves around the minimum. The scaling pa- 
rameters are chosen so that the elongations have approx- 
imately the same value in all directions. This amounts 
to an empirical form of preconditioning which makes the 
FES minimum nearly spherical in n dimensions and easy 
to fill with n-dimensional Gaussians. 

The metastep length 5a determines the efficiency of the 
method in a manner that is exponential in n and it should 
be made as large as possible. However, the dynamics in 
eqn. (|l|) is able to reconstruct details of the FES only on 
the length scale of 5a. Moreover, an over-large value of 
5a would cause the system to jump irregularly from one 
basin to the other. The value of 5a is chosen requiring the 
metatrajectory to remain localized in the FES minimum 
if the bias potential is not applied. With this choice of 
5a, a single step of the dynamics in eqn. (|l|) cannot lead 
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a system from a FES minimum to another, and the initial 
state at each new iteration can be generated from the last 
MD step of the previous iteration without creating major 
overlaps. Since the shake algorithm is used to impose a 
new set of cr*^^ values, large forces on the constraints are 
generated at the beginning of each microscopic dynamics, 
and the initial part of the trajectory has to be discarded 
for the calculation of F*~^^ . 

The height of the Gaussians W can be estimated as 
follows. We notice that we want to reach a situation 
in which the modified FES is flat. In such a case the 
forces coming from the FES and that coming from the 
Gaussian approximately balance each other. If we re- 
quire the maximum value of the Gaussian forces to be 
smaller than the typical FES force, we arrive at the rela- 

tion W/5a e~ 2 = al^Ff)^ with a < 1. In all cases so far 
studied a value of a close to .5 allows a fast escape from 
the local minima in the FES without significant loss of 
the underlying structure. 

In the future adaptive ways of determining all these 
parameters should be considered in order to adapt the 
procedure in an optimal manner to the local shape of the 
FES. 



II. RESULTS 

A. Dissociation of NaCl in water 

We first discuss an application of the method to the 
dissociation of a NaCl molecule in water. The most sta- 
ble state is the dissociated one, while the undissociated 
contact ion pair corresponds to a metastable local mini- 
mum. Here we shall study how the system escapes from 
this metastable state„We model the system using the 
AMBER95 force fieldEl, and solvate the NaCl complex 
in 215 TIP3P water molecules with periodic boundary 
conditions. The electrostatic interaction between classip 
cal atoms is taken into account by the P3M methodic. 
As collective coordinates, we use the distance VNa-ci be- 
tween Na and CI, and, in order to take into account the 
dynamics of the solvation shell during the dissociation, 
we also use the electric fields VNa and Vci on the Na 
and on the CI due to the water molecules within -^6.5 
A , of the ions. V^va and Vci depend significantly on the 
hydration pattern around the Na and on the CI ion. If, 
for example, a hydrogen bond to one of the two ions is 
formed or broken, the field on the ion changes by several 
kcal/mol. A dynamics of the form was performed on 
the system, with Sa = 0.25 and W — 0.3 kcal/mol. The 
scaling parameters As^ were 0.53 A, 32 and 32 kcal/mol 
for rpfa~ch ^Na and Vci respectively. The forces in eqn. 
(|l|) were evaluated by short MD runs with a time step of 
0.7 fs performed on 6 replicas of the system. The repli- 
cas were then equilibrated for 200 MD steps at 300 K, 
and the forces evaluated by averaging over the follow- 
ing 500 MD steps. Starting from a distance of 2.6 A, 
corresponding to a contact pair, the ions dissociate after 




FIG. 2: Free energy as a function of the collective coordi- 
nates TMa-ci, VNa and Vci for a Na-Cl complex in water. 
The two-dimensional free energy F2 shown here are obtained 
as F2 = — l//91og [J d2;exp(— /3F3)] , where (3 is the inverse 
temperatures, z is the coordinate not shown here, and the 
three-dimensional free energy F3 is estimated with eqn. 
The contours are plotted every 0.5 kcal/mol. The zero of the 
free energy corresponds to the metastable minimum of F in 
the contact ion pair (C in figure). 



120 coarse-grained iterations. The dynamics was then 
continued for another 350 steps imposing a maximum 
separation of the ion pair of 5 A, in order to study the 
structure of the free energy surface around the transition 
state (see Figure 2). Seven recrossings for the coarse- 
grained dynamics were observed. The transition state 
is located rj^a-ci ~ 4.02 A, VNa ~ —120 kcal/mol and 
Vci ^ 81.5 kcal/mol (see Figure 2). The overall topology 
of the FES confirms the importance of the solvent degrees 
of freedom for describing the reaction, since the transition 
region is in trapsverse orientation relative to the axis of 
the coordinate£il, and hence neither the distance nor the 
fields alone can provide an exhaustive description of the 
dissociation. The transition barrier (S Ih Figure 2), es- 
timated with thermodynamic integrationD along rNa-cii 
is of 3.3 kcal/mol. The one-dimensional free energy as a 
function of VNa-ci can also be obtained by integrating 
the three-dimensional FES (F3) with respect to VNa and 
Vci as Fi = -l//31og[/ dVNadVci exp{-m)]- This 
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leads to a barrier of 3.4 kcal/mol. Hence, the method is 
able to reproduce the one-dimensional free energy profile 
estimated with thermodynamic integration. 

We have verified that these collective coordinates iden- 
tify the transition state-and the reaction coordinate also 
in the dynamical sensellil. To this effect we prepared a 
set of replicas of the system with the critical values of 
the collective coordinates. We then followed the proce- 
dure described in RefEJ, assigning to the particles val- 
ues of the velocity chosen from a Maxwellian distribu- 
tion at 300 K, and following the trajectory backwards 
and forwards in time. By this procedure, we found that, 
within a time scale of 200 fs, the overwhelming majority 
of the replicas fell into one of the two basins of attrac- 
tion without any recrossing. Moreover, the time-reversed 
dynamics rleads systematically to the opposite basin of 
attractiorJlil. Finafe, if we prepare the systems with the 
SHAKE algorithm^ for a given value of the collective co- 
ordinates, and if the particle velocities are not reassigned 
before letting the system evolve freely, the time needed 
to go in either of the attraction basins is of the order of 
a few ps, and several recrossings are observed. This is 
due to the fact that the SHAKE algorithm imposes zero 
velocities in the direction of the constraint and it is only 
after thermalization, which takes place in a time scale of 
picoseconds, that the collective coordinates acquire suffi- 
cient velocity to evolve towards the FES minima. This is 
a strong confirmation that we have identified the correct 
transition state. 

In order to characterize the solvent structure during 
the dissociation, we have computed the Na- Water pair 
correlation function at the transition state S. The coor- 
dination number, obtained by integrating the pair cor- 
relation function up to a distance of 3 A, is 5.28. i-This 
value is consistent with the picture depicted in RefH, in 
which the sodium at the transition state is shown to be 
typically only five-fold coordinated. 



B. Isomerization of alanine dipeptide in water 

As a second application, we used our method to ex- 
plore the free energy surface of an alanine dipeptide _as 
a function of the backbone dihedral angles $ and ^'E3. 
ThCrSystem is described by the all-atom AMBER95 force 
fieldB and solvated in 287 TIP3P water molecules. The 
parameters of the coarse-grained dynamics of eqn. (|^) 
are 5a = 0.25 and W = 0.25 kcal/mol, while the scaling 
parameters Asi are 60° for both <I> and The other 
settings are the same as for the simulation of the Na- 
Cl complex. The FES of alanine dipept ide L n water has 
been studied in detail by other authorsE3tJ, who iden- 
tified several minima. Our dynamics (0) is able to cap- 
ture the overall features of the FES extremely quL 
Starting from a configuration close to the state Cregol 
the dynamics fills this attraction basin in approximately 
45 steps. The minimum is localized at $ ~ —80° and 
^ ~ 150°. The estimated free energy at the saddle point 
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FIG. 3: Free energy as a function of the backbone dihedral 
angles "I> and ^ for an alanine dipeptide in water after 40, 
100 and 200 steps of coarse-grained dynamics (|l|). The zero 
of the free energy corresponds to the the state Ur, located at 
"I> = — 63 and = —30. The contours are plotted every 0.2 
kcal/mol. The state Creq is located at $ = —70 and ^I* = 165 
and 1.3 kcal/mol. Finally, the transition state between Or and 
Creq is located at $ = —80 and ^ = 35 and 3.3 kcal/mol. If 
the dynamics is continued, less favourable areas of the phase 
space are explored (e.g. the region with $ > 0), but the 
free energy differences reported here are unchanged within 
the accuracy of the calculation. 



is 2 kcal/mol higher than the free energy in Creq- ^ 
compares to a value of 2.1 kcal/mol obtained in Ref. 
by thermodynamic integration. Our estimate is obtained 
with 45 steps of coarse-grained dynamics on 6 replicas, 
corresponding to a total of only ~ 120 ps of mycroscopic 
dynamics. After another 150 steps, the dynamics fills the 
attraction basin of the state ar, which corresponds to the 
equilibrium configuration of the system, and recrosses to 
Cjeq- Finally, the dynamics visits the $ > region of 
the configuration space after ^ 300 steps. In Figure 3 we 
plot the free energy as a function of $ and 5" after 40, 100 
and 200 steps of coarse-grained dynamics. Also small de- 
tails of the free energy landscape, like the location of the 
saddle points and the presence of secondary minima (e.g. 



5 



the state CsES), can be detected graphically after a few 
steps of coarse-grained dynamics. The one-dimensional 
free energy profile as a function of ^ obtained by analytic 
integration of the FES reproduces the profile obtained in 
Rcf.ca within 0.2 kcal/mol for all the values of 

Concerning the choice of the collective coordinates for 
this system, it should be remarked that the dihedral an- 
gles $ and ^' do not provide a copaplete description of the 
dialanine isomerization reactionca, and the real reaction 
coordinate should take into account also the solvent de- 
grees of freedom. Still, our results show that the method 
is able to reproduce the overall feature of the FES even 
if relevant degrees of freedom are not explicitly included 
in the collective coordinate space. This shows that the 
neglected degrees of freedom, although relevant for de- 
termining the reaction coordinate, are associated with 
relatively small barriers and are sampled efficiently dur- 
ing the microscopic dynamics for each value of $ and ^E", 
in spite of the relatively short runs. This is due to the 
use of the replicas and to the ability of the dynamics to 
retrace the same region of parameter space during the 
dynamics. 

III. DISCUSSION 

The method is based on ihe combination of the ideas 
of coarse-grained dynamical3 in the space defined by a 
few collective coordinatosj and on the introduction of a 
history-dependent biasaO. Constructing a dynamics on 
a free energy surface that depends on a few collective co- 
ordinates, we simplify enormously the complexity of the 
problem, which depends exponentially on the number of 
degrees of freedom. The FES is much smoother than 
the underlying potential energy surface and it is topo- 
logically simpler, with a greatly reduced number of local 
minima. The history-dependent bias prevents the system 
from vistiting regions it has already explored. This term 
is crucial for the efficiency of the method: for example, a 
Langevin dynamics in the space of the collective coordi- 
nates would provide scant information on the underlying 
FES and it would be slower in locating global minima. 
This has been checked by performing a Langevin dynam- 
ics for the dissociation of a NaCl molecule in water, with 
the same collective variables introduced in Section 2 and 
adding to eqn. (0) a Gaussi an noise at a temperature T 
that is gradually reducedcil. The number of evaluations 
of the forces required to converge to the minimum of the 
FES (i.e. the dissociated state) is well above 1000, if 
the temperature is decresed with a logarithic schedule, 
as required in order to epsure quasi-ergodicity in the col- 
lective coordinate spaccEj. A history dependent bias po- 
tential as defined in eqn. (^) but applied in a regular MD 
simulation without coarse graining the dynamics would 



be efficient in finding escapes, from the local minima, at 
it has been shown elsewhereEl, but it would not provide 
quantitative informations about the FES. In particular, 
equation (jsj) holds only because the dynamic is performed 
in the collective coordinates space using the derivatives 
of the free energy. Hence, both the coarse-grained dy- 
namics on the FES and the history-dependent bias are 
essential ingredients of the method, and only their com- 
bination allows an efficient and accurate determination 
of the FES. 

Another advantage of our search algorithm is that it 
can be easily parallelized over the P replicas, thus opti- 
mally exploiting present-day computer architectures and, 
moreover, does not require a very accurate determina- 
tion of the forces. In fact, in model calculations on an 
analytic PES, we have found that the algorithm is ef- 
fective even for noise levels as large as 20-30 % of the 
forces, since the height and the width of the gaussians 
is such that the coarse-grained dynamics explores a re- 
gion of space of the size of a Gaussian several times, and 
the system is only discouraged (and not prevented) from 
revisiting the same spot of configuration space. Hence, 
the errors in the forces tend to even out during the evo- 
lution. The dynamics (|]) is "self-healing", i.e. it is in 
principle capable to compensate the effect of complitely 
wrongly located gaussians or of a wrong additional bias 
potential. We can therefore use quite short runs and a 
small number of replicas to evaluate the forces. This is 
at variance with thermodynajpic integrationEEl and the 
blue moon ensemble methocoQ, where a very precise de- 
termination of the forces is needed and where treating 
multidimensional reaction coordinates is a daunting task. 
Other technique&y, like the umbrella sampling or bias po- 
tential methodsdi3'ElE3, allow an accurate estimate of the 
free energy surface in many dimensions. The drawback 
of these methods is that they are efficient only if a good 
approximation for the FES is known a priori, while this 
is not required in our method. 

A last advantage of the method is its ability to provide 
qualitative information on the free energy of a system in 
a very short time. For example, the overall topology 
of a free energy surface can be determined by very few 
coarse-grained dynamics steps using "large" Gaussians. 
Subsequently, our "qualitative" knowledge of the FES 
can be improved using smaller Gaussians, eventually re- 
ducing the dimensionality of the problem by exploiting 
the topological information obtained with the large Gaus- 
sians. 
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